{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "import sys\n",
    "import glob\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "import itk\n",
    "from itk import TubeTK as ttk\n",
    "from itkwidgets import view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CTP04.mha', 'CTP06.mha', 'CTP08.mha', 'CTP10.mha', 'CTP12.mha', 'CTP14.mha', 'CTP16.mha', 'CTP18.mha', 'CTP20.mha', 'CTP22.mha', 'CTP24.mha', 'CTP26.mha', 'CTP28.mha', 'CTP30.mha', 'CTP32.mha']\n"
     ]
    }
   ],
   "source": [
    "# NRRD Study Name\n",
    "studyname = '../Data/CTP'\n",
    "\n",
    "# NRRD Files\n",
    "directory = (studyname + '/')\n",
    "\n",
    "# Saved NRRD Files \n",
    "directory2 = (studyname + '-Reg/')\n",
    "if os.path.isdir(directory2) == False:\n",
    "    os.mkdir(directory2)\n",
    "    \n",
    "# Mask Creation and Location\n",
    "directory3 = (studyname + '-MinMax/')\n",
    "if os.path.isdir(directory3) == False:\n",
    "    os.mkdir(directory3)\n",
    "    \n",
    "pic_folder = os.listdir(directory)\n",
    "pic_folder = [pic_folder for pic_folder in pic_folder if \".mha\" in pic_folder]\n",
    "pic_folder.sort()\n",
    "print(pic_folder)\n",
    "num_images = len(pic_folder)\n",
    "\n",
    "im0Tmp = itk.imread(directory + pic_folder[int(num_images/2)], itk.F)\n",
    "\n",
    "resample = ttk.ResampleImage.New(Input=im0Tmp,MakeIsotropic=True)\n",
    "resample.Update()\n",
    "im0 = resample.GetOutput()\n",
    "immath = ttk.ImageMath.New(Input=im0)\n",
    "immath.Blur(1)\n",
    "im0Blur = immath.GetOutput()\n",
    "\n",
    "immath.Threshold(150, 800, 1, 0)\n",
    "immath.Dilate(10, 1, 0)\n",
    "mask0 = immath.GetOutputUChar()\n",
    "mask0Tmp = itk.GetArrayViewFromImage(mask0)\n",
    "mask0Tmp[0:4,:,:] = 0\n",
    "sizeZ = mask0Tmp.shape[0]\n",
    "mask0Tmp[sizeZ-4:sizeZ,:,:] = 0   #No need to update mask0 since mask0Tmp is a view of mask0 (shared memory)\n",
    "\n",
    "itk.imwrite(mask0, directory3 + 'mask.mha', compression=True)\n",
    "maskObj = itk.ImageMaskSpatialObject[3].New()\n",
    "maskObj.SetImage(mask0)\n",
    "maskObj.Update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "de9a85c79a9441d3a095648c155472fd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itk.itkImagePython.itkImageUC3; pr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(mask0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "***  7% : 64s : CTP04.mha  ***\n",
      "***  13% : 65s : CTP06.mha  ***\n",
      "***  20% : 63s : CTP08.mha  ***\n",
      "***  27% : 64s : CTP10.mha  ***\n",
      "***  33% : 60s : CTP12.mha  ***\n",
      "***  40% : 60s : CTP14.mha  ***\n",
      "***  47% : 62s : CTP16.mha  ***\n",
      "***  53% : 62s : CTP18.mha  ***\n",
      "***  60% : 61s : CTP20.mha  ***\n",
      "***  67% : 60s : CTP22.mha  ***\n",
      "***  73% : 61s : CTP24.mha  ***\n",
      "***  80% : 61s : CTP26.mha  ***\n",
      "***  87% : 59s : CTP28.mha  ***\n",
      "***  93% : 60s : CTP30.mha  ***\n",
      "***  100% : 61s : CTP32.mha  ***\n",
      "Done\n"
     ]
    }
   ],
   "source": [
    "Dimension = 3\n",
    "PixelType = itk.ctype('float')\n",
    "ImageType = itk.Image[PixelType, Dimension]\n",
    "\n",
    "imdatamax = itk.GetArrayFromImage(im0)\n",
    "imdatamin = imdatamax\n",
    "imdatamax2 = imdatamax\n",
    "imdatamin2 = imdatamax\n",
    "imdatamax3 = imdatamax\n",
    "imdatamin3 = imdatamax\n",
    "\n",
    "imFixedBlur = im0Blur\n",
    "\n",
    "for imNum in range(num_images):\n",
    "    start = time.time()\n",
    "    \n",
    "    imMoving = itk.imread( directory + pic_folder[imNum], itk.F )\n",
    "    \n",
    "    immath.SetInput(imMoving)\n",
    "    immath.Blur(1)\n",
    "    imMovingBlur = immath.GetOutput()\n",
    "    \n",
    "    imreg = ttk.RegisterImages[ImageType].New()\n",
    "    imreg.SetFixedImage(imFixedBlur)\n",
    "    imreg.SetMovingImage(imMovingBlur)\n",
    "    \n",
    "    imreg.SetRigidMaxIterations(3000)\n",
    "    imreg.SetRegistration(\"RIGID\")\n",
    "    imreg.SetExpectedOffsetMagnitude(20)\n",
    "    imreg.SetExpectedRotationMagnitude(0.3)\n",
    "    imreg.SetMetric(\"MEAN_SQUARED_ERROR_METRIC\")\n",
    "    \n",
    "    imreg.SetFixedImageMaskObject(maskObj)\n",
    "\n",
    "    #imreg.SetReportProgress(True)\n",
    "    \n",
    "    imreg.Update()\n",
    "    \n",
    "    tfm = imreg.GetCurrentMatrixTransform()\n",
    "    imMovingReg = imreg.ResampleImage(\"LINEAR_INTERPOLATION\", imMoving, tfm, -1024)\n",
    "    \n",
    "    itk.imwrite( imMovingReg, directory2 + pic_folder[imNum], compression=True )\n",
    "    \n",
    "    imdataTmp = itk.GetArrayFromImage(imMovingReg)\n",
    "    \n",
    "    imdatamax = np.maximum(imdatamax,imdataTmp)\n",
    "    imdatamin = np.minimum(imdatamin,imdataTmp)\n",
    "    imdataTmp[np.where(imdataTmp==imdatamax)] = 0\n",
    "    imdataTmp[np.where(imdataTmp==imdatamin)] = 0\n",
    "    imdatamax2 = np.maximum(imdatamax2,imdataTmp)\n",
    "    imdatamin2 = np.minimum(imdatamin2,imdataTmp)\n",
    "    imdataTmp[np.where(imdataTmp==imdatamax)] = 0\n",
    "    imdataTmp[np.where(imdataTmp==imdatamin)] = 0\n",
    "    imdatamax3 = np.maximum(imdatamax3,imdataTmp)\n",
    "    imdatamin3 = np.minimum(imdatamin3,imdataTmp)\n",
    "    \n",
    "    end = time.time()\n",
    "\n",
    "    percent = (imNum + 1) / num_images * 100\n",
    "    print('***  ' + str(round(percent)) + '% : ' + str(round(end-start)) + 's : ' + pic_folder[imNum] + '  ***')\n",
    "    \n",
    "print('Done')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = itk.GetImageFromArray(imdatamax3)\n",
    "out.CopyInformation(im0)\n",
    "itk.imwrite(out, (directory3 + 'max3.mha'), compression=True)\n",
    "\n",
    "out = itk.GetImageFromArray(imdatamin3)\n",
    "out.CopyInformation(im0)\n",
    "itk.imwrite(out, (directory3 + 'min3.mha'), compression=True)\n",
    "\n",
    "out = itk.GetImageFromArray(imdatamax3 - imdatamin3)\n",
    "out.CopyInformation(im0)\n",
    "itk.imwrite(out, (directory3 + 'diff3.mha'), compression=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = itk.GetImageFromArray(imdatamax)\n",
    "out.CopyInformation(im0)\n",
    "itk.imwrite(out, (directory3 + 'max.mha'), compression=True)\n",
    "\n",
    "out = itk.GetImageFromArray(imdatamin)\n",
    "out.CopyInformation(im0)\n",
    "itk.imwrite(out, (directory3 + 'min.mha'), compression=True)\n",
    "\n",
    "out = itk.GetImageFromArray(imdatamax - imdatamin)\n",
    "out.CopyInformation(im0)\n",
    "itk.imwrite(out, (directory3 + 'diff.mha'), compression=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "460952a317fa4680865b999f6140752c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itk.itkImagePython.itkImageF3; pro…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "out = itk.GetImageFromArray(imdatamax-imdatamin)\n",
    "out.CopyInformation(im0)\n",
    "view(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
